Analysis of the quantum-classical Liouville equation in the mapping basis 
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The quantum-classical Liouville equation provides a description of the dynamics of a quantum 
subsystem coupled to a classical environment. Representing this equation in the mapping basis 
leads to a continuous description of discrete quantum states of the subsystem and may provide 
an alternate route to the construction of simulation schemes. In the mapping basis the quantum- 
classical Liouville equation consists of a Poisson bracket contribution and a more complex term. By 
transforming the evolution equation, term-by-term, back to the subsystem basis, the complex term 
(excess coupling term) is identified as being due to a fraction of the back reaction of the quantum 
subsystem on its environment. A simple approximation to quantum-classical Liouville dynamics in 
the mapping basis is obtained by retaining only the Poisson bracket contribution. This approximate 
mapping form of the quantum-classical Liouville equation can be simulated easily by Newtonian 
trajectories. We provide an analysis of the effects of neglecting the presence of the excess coupling 
term on the expectation values of various types of observables. Calculations are carried out on 
nonadiabatic population and quantum coherence dynamics for curve crossing models. For these 
observables, the effects of the excess coupling term enter indirectly in the computation and good 
estimates are obtained with the simplified propagation. 
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I. INTRODUCTION 

The construction of algorithms for simulating the 
quantum dynamics of an arbitrary many-body system 
is a long standing problem. Since the computational 
cost of a quantum simulation scales exponentially with 
the system size, a full quantum calculation is impracti- 
cal except for small systems. This fact has prompted 
the construction of a variety of approximate methods for 
the simulation of quantum dynamics. One class of ap- 
proximate schemes singles out a portion of the system 
for a full quantum treatment while its environment — the 
rest of the system — is treated classically. A number of 
such mixed quantum-classical schemes have been pro- 
posed and references to this literature can be found in 
reviews on this topioi~— . 

We focus on one scheme of this type based on the 
quantum-classical Liouville equation (QCLE), 



d_ 

di 



p w (X,t) 



[H w ,pw(t)] 



(1) 



({H w ,Pw(t)} - {pw(t),Hw}), 



where [•, •] is the commutator and {•, •} is the Poisson 
bracket. Here X = (X x , X 2 , . . . , X Na ) = (R,P) = 
(Ri, R%, ■ ■ ■ j Rn s , Pi, Pi, ■ ■ ■ , -P/v e ) are the positions and 
momenta of the 2N e environmental degrees of freedom 
and the index W stands for a partial Wigner transform 6,7 
over the environmental degrees of freedom, so that the 
density matrix, pw(X), and Hamiltonian, H\y(X), are 
operators in the Hilbert space of the subsystem while 
functions of the phase space variables for the bath de- 
grees of freedom. (The dependence on the phase space 
coordinates will be omitted when confusion is unlikely to 
arise.) The nature of QCL dynamics and the statistical 
mechanics of systems following this dynamics have been 



described^. For a review with references to the literature 
on this topic, see Ref. ||. 

While the QCLE has a number of attractive features 
and has been shown to provide an accurate description 
of the dynamics in many instances*^—, it is difficult to 
solve. A direct numerical integration of Eq. (|TJ) requires 
very fine spatial grids and can be used for the study of 
small systemaii. The QCLE can be cast in any basis 
which spans the Hilbert space of the quantum subsys- 
tem^. When written in the quantum subsystem basis, 
the QCLE can be solved using a trajectory-based algo- 
rithm where trajectories are not independent of one an- 
other^i^. Representation of QCLE in the adiabatic basis 
gives a more intuitive picture in terms of classical tra- 
jectories moving on single or mean Born-Oppenheimer 
surfaces^— . Also, diagonalizing the Hellman-Feynman 
force derived from the Born-Oppenheimer potential leads 
to the force basis which yields yet a different route to 
simulation of the QCLE^ZL. 

Another representation of QCL dynamics is in terms 
of the mapping basi o 28 i 29 . This basis, which provides a 
description of a discrete quantum system in continuous 
variables, has been used in a number of different appli- 
cations to quantum dynamical problems^ - — . When the 
QCLE is expressed in this basis one obtains an evolu- 
tion equation where the evolution operator consists of 
a Poisson bracket contribution which can be solved by 
characteristics, and a more complex term that involves 
specific correlations of the quantum subsystem and clas- 
sical environment. Thus far, simulations of the QCLE 
in the mapping basis have neglected this more complex 
contribution. In this paper we provide an analysis of the 
mapping form of the QCLE that elucidates the nature 
of this neglected term. By transforming the QCLE in 
the mapping basis back to the subsystem basis, where 
an interpretation of the physical meaning of the different 
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terms is easier, we show that the neglected term (called 
the excess coupling term) accounts for a fraction of back 
reaction of the quantum subsystem on its environment. 
We also describe how the contribution of the excess cou- 
pling term to the average values of certain observables 
can be estimated numerically. 

The outline of the paper is as follows: In Sec. |TT] we 
summarize the equations to be analyzed and introduce 
some key definitions, while in Sec. Mil we perform the 
term-by-term analysis of the mapping equation based on 
the back transformation to the subsystem basis. Sec- 
tion IIVI presents an analysis of the approximate form 
of the mapping QCLE that includes only the Poisson 
bracket term, called the Poisson bracket mapping equa- 
tion (PBME) . This section also describes how neglecting 
the excess coupling term affects average values of differ- 
ent types of observables. In Sec. E] the performance of 
the PBME is tested on two often-studied curve crossing 
models, and the effects of the excess coupling term in 
the evolution operator are discussed. The last section 
contains the conclusions of our investigation. 



II. SUBSYSTEM AND MAPPING BASIS 
REPRESENTATIONS 

For a system with partially Wigner transformed Hamil- 
tonian H w = £L + £- + V s (r) + V e (R) + V C (R, f), where 

V e , V s and V c are, respectively, the environment, subsys- 
tem and coupling potentials, R and P are again the N e - 
dimensional coordinates and momenta of particles of the 
environment with mass M , and f — (fi, . . . , tn s ) and 
p = {pi, p2, ■ . ■ , Pn s ) are the N s coordinate and momen- 
tum operators of the particles of the quantum subsystem 
with mass m, the QCLE in the subsystem basis takes the 
form—, 

^ P ^'(X,t) = -iLJ aa ,p<^'(t) 

+\{pw" wf a ' -vr" Pw a ' (t)) 

fdV e d P d \ 
+ \M'dP~M'dR) P ^ {t) 
1/ dVr" dp^'jt) dp<$"{t) dVf a ' \ 
2\ dR dP dP OR J 

= -iC aa ,^ P ^'(X,t). (2) 

Here and in the following we use the Einstein sum- 
mation convention where repeated indices are summed. 
The subsystem basis is defined by the eigenvalue prob- 
lem, h s \a) = e a \a), where h s = ^ + V s . In Eq. ©, 

u aa > = (e a - e a ,)/n and V c aa ' = (a|F c |a'}- The last 
line in Eq. ^ defines >C QCt ', MM ', the QCL operator in the 
subsystem basis. 

Starting from this subsystem representation of the 
QCLE, we may transform it to the mapping basis in 



the following way. Suppose that there are N subsys- 
tem quantum states |A), A = 1, . . . , N. These subsystem 
states may be mapped onto harmonic oscillator states via 
the transformation, |A) — > \m\) = | Oi , • - • , 1a, • • • ,0jv), 
where 



(q\mx) = (qi,q 2 , 
= <M<7i) 



■ ! Qn\0i, ■ ■ ■ , 1a, • ' ' i Ojv) 
•<A)(?a-i)<M?a)---<Ao(?jv), (3) 



with <pQ and <p\, respectively, being the ground and the 
first excited state wave functions of a harmonic oscillator. 
The creation and annihilation operators on the mapping 
states, cl\ = (q\ + ip\)/V2h and a\ = (q\ — ip\)/\/2h, 
are used to define a mapping of operators as^ 



A = A XX |A)(A'| -¥ A Tl 



A xx 'a\a w . 



(4) 



For example, the partially Wigner transformed mapping 
form of the system Hamiltonian is 

P 2 1 

— + V e (R) + —h xx (R)(Mx' +P\P\> - 



2M 



V e (R) + h m (R). 



(5) 



Here we used the expression for the creation and annihi- 
lation operators in terms of the coordinates and momenta 
of the harmonic oscillators, together with the appropri- 
ate commutation relationship. Also, h xx (R) — (X\h\X') 
with h = h s + V C (R, f) so that, more explicitly, 



h xx '(R) = e x 5 xy +V XX \R). 



XX' i 



(6) 



In writing Eq. ([5]), we assumed that h xx = h x x . Follow- 
ing this prescription, after a further Wigner transforma- 
tion^ over the mapping variables, 



p m {X, x) 



(2-Kh) 



N 



dz ^{r-J^WIr+S 



(7) 

where x = (x±, x%, . . . , xn) = (r,p) — 
(ri, r 2 , . . . , r N ,pi,p 2} . . . ,pn), the QCLE takes the 
form^9 

h xx ' 



i-p m (X,x,t) = - — 



fdHra 

V dR 



dh 



_d_ 



px, -k- rx, ik) Pm{t) (8) 

P d 
M ' dRJ 



)Pm(t) 



d 2 



d 2 



dR \drxdrx 1 dpxdpx 



dp m (t) 

dP 



Note that the Wigner transform variables r, p and z have 
the same dimension as the number of quantum states, N . 
The Wigner transform of H m over mapping variables is 



H m (X,r,p) = J dz jP^( r -~\H n 



^ + V e (R) + —h xx '{R)(rxrx> +P\P\> - K5 X x- 



p2 

2M 
2M 



V e (R) + h m (R). 



(9) 
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We can also write the mapping form of the QCLE as 
d 



Of 



p m (X,X,t) = {H m ,p m (t)} X ,x 



dh 



A A' 



d 2 



d 2 



dR \dr\dr\i 



dp x dp x 



(10) 

dp m {t) 

dP 



In Eq. (fTOj) . {H m , p m }x,x denotes a Poisson bracket with 
respect to both the bath X and mapping x variables. Ne- 
glecting the last term of the dynamics in Eq. (|10p yields 
a Hamiltonian system of equations which can be easily 
solved using Newtonian trajectories. The last term has 
a complex form involving both bath and mapping dif- 
ferential operators which make its interpretation difficult 
and precludes the implementation of simple algorithms 
for the simulation of the full mapping form of the QCLE. 



III. TRANSFORMING MAPPING DYNAMICS 
BACK TO THE SUBSYSTEM BASIS 

In order to understand the nature of the last term in 
the mapping QCLE, we shall transform each term in this 
equation back to the quantum subsystem basis. Natu- 
rally, combining all the back-transformed terms, the orig- 
inal subsystem basis equation (Eq. ([2])) is recovered; how- 
ever, the term-by-term transformation highlights how 
each contribution in the mapping representation of the 
QCL operator is associated with a specific contribution 
in the subsystem basis form of the operator. This proce- 
dure leads to a simple physical interpretation of the last 
term in Eq. (|TU)) . 

We begin by recalling the expression for the Wigner 
transformed density matrix in the mapping basis. Us- 
ing the definition of a mapping operator in Eq. Q and 
the expression for the Wigner transform of the density 
operator in Eq. (J7J), one has 



Pm(X,x) = 



1 



{2irh) 



N 



dze^' h (r--\p m (X)\r+-) 



1 



r2d/; , J dze^/ h p%\x){r~-\a\ay\r+-) 
p^'(X)c xx ,(x), (11) 



where 



CAA' (x) 



1 



2h(2irh) N 

x r x r x , + i(r\p\> - r x >p\) + p\P\> - %5\\> 



(12) 



Details of the calculation of caa' (%) are given in the Ap- 
pendix. 

Next, given p m (X,x) we consider how to recover the 
expression for pfy (X). The mapping relationship in 
Eq. (j4|) ensures that the matrix elements of an opera- 
tor are the same in the subsystem and mapping bases, 
thus, 



P^'(X) = (X\ Pw \X') = (m x \p m \mx< 



(13) 



Inserting resolutions of the identity in the coordinate rep- 
resentation, the last term in the equality above can be 
written as 

(m\\p m \mx>) = / dydy' (m x \y)(y\p m \y')(y'\m x ,) (14) 



drdz (m x \r ~^){r- ^\p m \r + -}{r + -\m x >) 



Mean r = (y + y')/2 and difference z = y' — y coordi- 
nates were introduced in the last line to pave the way for 
the introduction of the Wigner transform of the density 
operator in the mapping coordinates. The off-diagonal 
element of the operator in the equation above can in fact 
be expressed as (r — ||/5 m |r+ |) = J dpe~ lp ' z / h p m (X,x) 
(the inverse transform of the expression in Eq. ([?]))■ Com- 
bining this expression with the identity in Eq. (|13[) . we 
get 



Pw'{X) 



(15) 



drdzdpe lp - z/h p m (X,x)(m x \r- -)(r+-\my) 



In the Appendix we show that the integral over z in the 
equation above can be performed analytically to obtain 



(${X)= I dx p m (X,x)g xx ,(x), 



(16) 



where 



2 N+1 



-x 2 /h 



(17) 



r X r x > - i(r x p x > - r x ,p x ) + p x p x > - -5 XX , 



(In the notation of this paper x 2 = x ■ x — r 2 + p 2 — 
^2 x (r x +p\)- For clarity, in some places we shall use the 
more expanded forms of this condensed notation.) Rela- 
tionships analogous to those in Eqs. (ITT1) and (|T6| hold 
for the representation of any operator in the two bases, 
except for the fact that the multiplicative prefactors in 
Eqs. (fT2")) and (fTT|) are interchanged. 

These results can be used to establish the relation- 
ship between the mapping and the subsystem bases. To 
set the stage for our analysis, we multiply each term of 
Eq. (p~0|) by g aa '{x) and integrate over the mapping co- 
ordinates to obtain 







Q~ t Pw C*'*) = / dx g aa ,{x){H mi p m {t)} x ,x 



dx g aa , (x) 



dh 



A A' 



d 2 



OR \dr x dr X i dp x dp x 



d_ 



XPm{t). 



(18) 



In the expression above, the integrand contains 
p m (X, x, t). We can use Eq. (fTTj) to express this quantity 
in terms of the density matrix in the subsystem basis, 
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Pw (*)' t0 § et 
d 

-Q t Pw PM) = I dx 9aa'{x){Hm,Pw (t)Cw'( X )}x,x 



This integral too can be performed after substitution of 
the g aa > and c MM < functions to get 



dx g aa > (x) 



xp$(t)c^(x). 



dh 



A A' 



(- 



if 



d' 2 



dR \dr\dr\i dp\dp\> 



d_ 

dP 



(19) 



In order to proceed with our analysis, we consider sepa- 
rately the Poisson bracket and second complex terms on 
the right hand side. 



A. Transformation of the Poisson bracket term 

The hrst contribution in the Poisson bracket, 
{H m , p m (t)}x,x, arising from derivatives with respect 

to the mapping variables is (p\> afj - »> g§H Pm 

(see Eq. ©). Its back transformation is given by 

h xx \ ( d d 



PBl 



drdp 



PX' 



fifi 
Pw C P 



h ' dr x dp X/ 

(20) 

Direct substitution of the expressions for g aa i and / w / 
results in a sum of integrals that can be evaluated ana- 
lytically (they all are in the form of a polynomial times 
a Gaussian function) to obtain 

PBl = -iCj aa >Pw + l j t {pw"yf a - V c aa "pw a ) ( 21 ) 

after substitution of h aa from Eq. ©. This contribution 
represents the evolution of the quantum subsystem and 
includes the influence of the environmental degrees of 
freedom through the coupling potential V c . 

Next, we consider the back transformation of the sec- 
ond contribution in Eq. (jSJ coming from derivatives with 
respect to environmental coordinates, 
mi)Pm(t)- Given the form of H m in Eq. ^ 
write this contribution as 



d 

dP 



P_ 

M 

one can 



dH m dpr, 



dR OP 



P 

M 



dPm 

OR 



dPm 

dP 

dp„ 



dVe 
OR 

dh m 
dR dP 



P 
M 



dp m 

dR 
(22) 



Because dV e /dR = ~F e (R) and P/M are independent 
of the mapping variables, the back transformation of the 
first two terms in Eq. (|22l) is simple: 

/fdV d P d \ i 
drd ^A-dk-d-p~M-^ c ^ 



d_ 
~dP 



P 

M 



d_ 

dR, 



)pw 



(23) 



The last term in Eq. (|2"2"|) , arising from ^frr, gives 



PB3 = J drdp g aa ,-^-(r x r x , 



+ P\Py - hS\\>) 



dh xx dp^_ 
~dR d~P 



(24) 



PBS 



= if* 



l/dV c aa " d P ^ a ' dp^" dV c a " a 



2\ dR 



J w 
dP 



w 
dP 



dR 



l^(9Vo.dPw\ 5 , 
4 \ dR dPJ aa ' 



(25) 



where Tr' indicates a trace over the quantum subsystem 
states. 

Combining the results above, we obtain 



dx Q&a,' (^) { -^m ? Pm (^) } y- x i^aa' P\y (^) 



1 r V j( dV o dp W \ x 

"4 \dR' dP J aa ' 



(26) 



Thus, we find that the back transformation of the Pois- 
son bracket term yields the QCL operator in the subsys- 
tem basis defined in Eq. ([2]), plus an extra contribution 

Z^ 1 ' (lm ' ^§PJ^ aa '' K is this term that is responsi- 
ble for any errors in the approximate simulations of the 
QCLE that include only the Poisson bracket term. 

An understanding of the physical meaning of this extra 
contribution can be obtained from the following consid- 
erations: First, we observe that the trace of the partially 
Wigner transformed density matrix over the quantum 
subsystem states is just the phase space density of the 
environment, pb(X,t) = Tr' pw(X, t). Taking the trace 
of Eq © yields, 



or, equivalently, 

T ,(dVc dp w (t) 
\dR' dP 



d P d \ 
Ml' %^)' ™ 



d_ 

dP 



Fe(R) 

Pb(t), 



P 

M 



d_ 

dR 



Pb(t) 
(28) 



where iL e is the classical Liouville operator for the envi- 
ronment in isolation from the quantum subsystem. Con- 
sequently, the trace term in Eq. (|28|) can be interpreted 
as the time variation of the phase space density of the en- 
vironment along the flow lines generated by the classical 
evolution of the environment in isolation from the quan- 
tum subsystem. In a system where there is no coupling 
between the environment and the quantum subsystem, 
this term is zero by Liouville's theorem. As a result, 
a non-zero value of the trace contribution can be inter- 
preted as an effect arising from the back reaction of the 
quantum subsystem on its environment. 



B. Excess coupling term 

To complete the back transformation of the mapping 
QCLE to the subsystem basis, we must back transform 
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the complex last term in Eq. ([TO 



- / drdp g a 



dh 



XX' 



if 



dR dr\dr x 



dpxdpx 



_d_ 

dP 



x Pw c w' 



dh 



xx' 



dp 



xx 
w 



dh 



XX' 



dp 



A' A 
W 



dR 

,dV c 



4 TT '(- d R 



dP 

dpw - 
dP ' 



dR 



dP 



(29) 



Since this result is again proportional to Tr'^^^- • 
dpwjt) ^ ^ we ca j[ this contribution a excess coupling term. 

Inserting these results into Eq. (TT9"|) . the QCLE in the 
subsystem basis (Eq. ©) is obtained as expected. How- 
ever, from this derivation, we learn that the excess cou- 
pling term is equal and opposite in sign to the last term 
in Eq. (|26|) . As a result, this contribution exactly can- 
cels the analogous term in the back-transformed Poisson 
bracket expression to yield the desired result. 



IV. APPROXIMATE MAPPING QCLE 

The results of the previous section can be used to as- 
sess the utility of the approximate mapping QCLE where 
only the Poisson bracket contribution is retained in the 
Liouville operator. This Poisson bracket mapping equa- 
tion (PBME) is 



d_ 

di 



Pm(t) ~ {H m , p m (t)}x,x 



+ ■ 



n 

dH m 
~dR 



d 



' py ik ry Wx ]l 



dpm{t) P dp m (t) 



dP 



M dR 



(30) 



It can be solved by characteristics^ and the resulting set 
of ordinary differential equations is^ 



dr x {t) 

dt 
dpxjt) 

dt 
dR(t) 

dt 



1 



h xx '(R(t))px<(t), 



1 

m 

M 



h xx '(R(t))r y (t), 

dP(t) _ 
dt 



dH m 
~dR(ty 



(31) 



Consequently, the simulation of the dynamics described 
by Eq. ([30]) is an easy task4i^ 

From the results in Sec. [TIT] it follows that Eq. ([3"0")) 
is equivalent to the following equation in the subsystem 
basis: 



d 

—Pw {X,t) 



iC aa ' ,fj,fj,' Pyy (t) 

I ,(dV c dpw(t) \ 

-i Tr {-dR--dp-) 5aa '- (32) 



Because of the presence of the second term on the right 
side, this form of the equation shows that the back reac- 
tion of the quantum subsystem on the dynamics of the 
environmental degrees of freedom is incorrectly described 
and the solution by characteristics is only an approxima- 
tion to full QCL dynamics. Two features should be kept 
in mind when considering the error incurred in the use of 
the PBME. First, the effects of the environment on the 
quantum subsystem are fully accounted for in the Pois- 
son bracket, and the effects of the quantum subsystem on 
the environment are also included in this term through 
the first term in PB3 in Eq. ([2"5"f. Second, note the factor 
of 1/4 in Eq. ([321) : the error involves only a portion of 
the back reaction of the quantum subsystem on the evo- 
lution of the environmental density. To the extent that 
the effects of the excess coupling term are small, simula- 
tions employing the PBME for the evolution will provide 
an accurate description of the dynamics. 

The equivalence of Eqs. f[32]) and (f30|) suggests a means 
to gauge the importance of the excess coupling term. 
Rather than focussing on the excess coupling term it- 
self, it is more convenient to consider its effect on the 
estimates of expectation values of observables. The ex- 
pectation value of an arbitrary operator B\y{X) can be 
written as^ 

W) = J dX B^'(X)p^ x (X,t) 



= I dX (X,t)p^ x (X) 

dXdx B m (x, X, t)p m (x, X), (33) 

In the first line of Eq. ([3"3")l the expectation value is com- 
puted in the subsystem basis, in the second line the time 
evolution is moved from the density matrix to the oper- 
ator, while the last line gives the expectation value com- 
puted in the mapping basis. The quantities entering in 
the mapping form of the expectation value are the time 
dependent observable, 

B m (x,X,t) = {2h)~ 1 B xx ' (X(t)){r x {t)r x ,{t) (34) 
+px(t)px'(t) +i(rx(t)p x >(t) -p x (t)r x >(t)) - hS X x 



and the initial density p rn (x, X) = J dx'f(x, x')p m (x\ X) 
with 



f(x,x') 



1 



(27^)" 

z 1 



dzdz' (mx\r )(r H — \my) 



x{mx>\r' - ^}{r' + -\mx)e'' (p - z+p '- z "> /h . (35) 

The integrals involved in the definition of the function 
f(x,x') in Eq. ([351 can be performed analytically. For a 
two-level system the result is 



f(x,x') 



ir 2 h 4 



2(r • r' +p ■ p' f + 2(r -p' - r' ■ pf 



-h(x 2 +x' 2 ) + h 2 }e-^ 2+x ' 2 ^ h 
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where r ■ r' = rir[ + T2r' 2 , with similar expressions for 
other scalar products. Thus, these formulas allow one to 
use either the mapping or subsystem bases to compute 
the expectation values of any observable. 

While Eq. ([55]). along with Eqs. (JM]) and ([SI]), provide 
a simple means to compute the expectation value of an 
operator using the PBME, it is convenient to consider 
the computation of the general operator B\y(X) using 
the PBME written in the subsystem basis in order to 
gain further insight into the effects of the excess coupling 
term on such average values. Taking the time derivative 
of Eq. fl33j and using Eq. ([32]) we find 



-% J dX B^ a (X)£ aa ,^p^ (t) (36) 



dB(t) 
dt 



= -i J dX B^ a {X)C aal ,^p^ (t) 

1 [ dB&{X) dV c ^' VA 
~Zj dX ^P dR~ Pw W<W ' 

where, as usual, the Einstein summation convention is 
used. In the last line of this equation an integration by 
parts with respect to P was performed so that the mo- 
mentum derivative of the density matrix no longer ap- 
pears. Therefore, the effect of the excess coupling term 
on the evolution of the expectation value of Bw(X) can 
be estimated by computing the average values on the 
right side of this equation and determining their relative 
magnitudes. 

If the initial value of the operator Bw(X) is a function 
only of the configuration space coordinates, or indepen- 
dent of environmental phase space coordinates, then the 
last excess coupling term in Eq. (|36[) will vanish because 
the momentum derivative is zero, and for this case we 
have 



of these observables are independent of the phase space 
coordinates. 

If instead the observable is a function G(P) only of 
the environmental momenta, then the excess coupling 
term enters the computation of the time rate of change of 
its expectation value directly. Consider the expectation 
G{P)8 aa ,, 



value of B^ a (X) 



G W=E / dX G(P)p^(X,t). 



(39) 



If we compute the time derivative of this quantity using 
Eq. ([36]) we obtain, 



dG(t) 
dt 



dX^p--F aa 'pi a {X,t) (40) 



where F aa = -d(V e {R) + V c aa {R))/dR is the total force 
on the environment and F C QQ ' = -dV^ a ' (R)/dR is the 
force contribution arising from the quantum subsystem- 
environment coupling. The first contribution on the 
right side of Eq. (|4T)]) is obtained from the explicit 

computation of — i j dX G(P)6 a i a £ aa > . w ,< p 1 ^ (t), while 
the second contribution comes from the evaluation of 
\jdX G(P)S a , a Tr'(§? ■ ^iy aa ,. These formu- 

las provide an indication of how the time evolution of 
the environmental momenta are influenced by the excess 
coupling term. In particular, if the coupling is weak or if 
the environment is large and only a portion of the envi- 
ronment is directly coupled to the quantum subsystem, 
the total force on the environment will dominate the cou- 
pling force and this effect will be small. 



CURVE CROSSING MODELS 



dB(t) 
dt 



= -i dX B 



Pw (*)■ 



(37) 



For such observables the excess coupling term will enter 
the computation of the average value indirectly through 
the density matrix on the right side. This can be seen by 
writing the formal solution of Eq. (]3"2"]l as 



Pw'(t) = (e-^ f ) ftf(0) 

\ / act 



(38) 



4 OR dP 



and inserting the result on the right of Eq. ([57]) . The eval- 
uation of this higher order correction is difficult because 
of the convolution in the above equation. In particular 
these considerations apply to the computation of impor- 
tant observables such as the quantum subsystem popu- 
lations and quantum coherence, since the initial values 



In this section we consider the simulation of quantum- 
classical Liouville dynamics in the mapping basis for two 
curve crossing models studied earlier by Tully*£. In 
particular we carry out simulations using the PBME 
and compare the results with numerically exact quan- 
tum dynamics using a discrete Fourier transformation 
method^ifor these systems. Previous studies of the spin- 
boson model showed that the PBME yields results which 
are in excellent agreement with the numerically exact 
quantum simulations for the entire range of system pa- 
rameters that were studied^. In the spin-boson model, 
the coupling between the environment and the quantum 
subsystem is linear in the environmental coordinates and 
is non-zero for all values of this quantity. The curve 
crossing models studied here provide a further test of 
the PBME for a situation where the coupling among the 
different degrees of freedom is nonlinear in the environ- 
mental coordinates and is localized in the configuration 
space of the environment. 
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A: Simple avoided crossing 

The simple two-state avoided crossing is defined by the 
following Hamiltonian matrix in the diabatic basis^: 



A[l 



_ P -B\R\]R 
\\R\ 



Ce 



-DR* 



Ce- DR2 



(41) 



The corresponding total Hamiltonian in the mapping ba- 
sis is 



p2 


i 


2M 


" 2h 


C 


-DR- 


+r 





R 



4i --"W 



(r x r 2 + pip 2 )- 



2 2 2\ 

Pi~r 2 -p 2 ) 
(42) 



The diagonal elements of the Hamiltonian in the dia- 
batic basis and the adiabatic energies for this model are 
sketched in Fig. [T] Our simulations are carried out using 
atomic units (a.u.) and the parameters in these units are 
taken to be A = 0.01, B = 1.6, C = 0.005 and D = 1.0. 

We assume that the initial density matrix is uncor- 
related and can be factored in to a product of subsys- 
tem and environmental contributions, i.e., p m (X,x) — 
p s (x)p e {X), where p s {x) is the subsystem density matrix 
in the mapping basis and p e (X) is the distribution func- 
tion for the environment. The density p e {X) is chosen 
to be the Wigner transform of a Gaussian wave packet 
centered at i?o with momentum Pq: 



Pe(X) 



1 



,-{R-R ) 2 /<J 2 



(43) 



The mass M of the environmental degree of freedom is 
taken to be 2000 atomic mass units, while i?o = —3.8 
a.u. and a = 1. 

We assume the subsystem is initially in state 1 , so that 
in the mapping basis 



Ps(x) 



Pi 



(44) 



We consider the time evolution of the subsystem popu- 
lations, p" a (t), a = 1,2, and coherence as measured by 
the off-diagonal element of the subsystem density matrix, 
pl 2 {t). These quantities are conveniently computed us- 
ing the formulas in Eq. (|33j) for a general operator by 
selecting B$' = S\ a S\* a and B$' = S\i5\> 2 , respec- 
tively. The asymptotic values of the populations in states 
1 and 2 after the system has passed through the inter- 
action region are shown in Fig. [1] as a function of the 
initial momentum of the wave packet. We see that the 
PBME and exact quantum results are in excellent agree- 
ment for values of Pq above 10. Under this threshold, 
nuclear quantum effects play an important role so the 
most likely cause for the discrepancies in the figure is the 
breakdown of the classical dynamics approximation for 
the nuclei (note that usually, mixed quantum-classical 
non-adiabatic tests do not explore the range below this 
value of the momentum for this model). 





FIG. 1: (top) Diagonal elements of the Hamiltonian in the di- 
abatic basis (dashed lines) and adiabatic (solid lines) energies 
for the simple avoided crossing model, (bottom) Asymptotic 
populations of state 1 (solid squares with error bars) and state 
2 (solid circles with error bars) as a function of the initial mo- 
mentum of the wave packet, Po, for the simple avoided cross- 
ing. The corresponding numerically exact quantum results 
are indicated by open squares for state 1 and open circles for 
state 2. The exact results are connected with lines as a guide 
for the eye. 



The time evolution of the real and imaginary parts 
of pl 2 (t) are displayed in Fig. [5] These results for the 
quantum coherence are also in good agreement with the 
numerically exact full quantum simulations. The com- 
parisons shown in these figures test two effects at the 
same time: the validity of the QCLE and its approxima- 
tion by the PBME. Recall that from the considerations 
in Sec. |IV] the effects of the excess coupling term enter 
the populations and quantum coherence only as higher 
order contributions. 

For the simple avoided crossing (and the dual avoided 
crossing considered below), V e (R) = so that F aa = 
F" a and the two integrals on the right of Eq. (|40|) are 
equal, and the two contributions differ only in their pref- 
actors. For this two-level model N = 2 and the PBME 
prediction for the time rate of change of G(t) is 3/2 that 
of QCL dynamics. This is a severe test of the approx- 
imation since V e (R) — for this model. Thus, there is 
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FIG. 2: The real (thick lines) and imaginary (thin lines) parts 
of pl 2 (t) versus time for the simple avoided crossing model. 
The dashed lines denote the full quantum results while the 
solid lines are the results of computations using the PBME. 
(top) P = 50, (bottom) P = 10. 



FIG. 3: (top) Diabatic diagonal elements of the Hamiltonian 
(dashed lines) and adiabatic energies (solid lines) for the dual 
avoided crossing model, (bottom) Populations of the diabatic 
states with same symbols as those in Fig. [T] 



a substantial influence on the environmental momenta 
but, as noted above, these enter only in higher order 
corrections to the quantum subsystem populations and 
coherences, contributing to the good agreement with the 
exact results seen in the figures. 

B: Dual avoided crossing 

The dual avoided crossing model provides a further test 
of the theory. The Hamiltonian matrix in the diabatic 
basis takes the form^*, 

resulting in the mapping Hamiltonian, 

Hm = St + "k ( - Ae ~ BIi ' 1 + E °)irt+Pl-1) 

+ ^-e- DR \r 1 r 2 + PlP2 ). (46) 

The system parameters were taken to be A = 0.10, 
B = 0.28, C = 0.015, E Q = 0.05 and D = 0.06, again in 



atomic units. The initial conditions have the same form 
as for the simple avoided crossing, with Rq = —10. The 
diagonal matrix elements in the diabatic representation 
and the adiabatic energies for this model are plotted in 
Fig. [3] (top) and the populations as function of Pq in the 
bottom panel of this figure. One can see that the agree- 
ment with exact quantum dynamics is very good over 
the entire range of Po, although there are small discrep- 
ancies in the magnitudes of maxima and minima in the 
population curves. 

The time evolution of the real and imaginary parts 
of pl 2 {t) for the dual avoided crossing model are shown 
in Fig. |4] The results for the quantum coherence are 
again in good agreement with the numerically exact full 
quantum simulations for this model. 



VI. CONCLUSION 

By transforming the mapping form of the QCLE term- 
by-term back to the subsystem basis we were able to 
identify the nature of the complex term in this equation 
that makes its simulation difficult. This observation then 
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FIG. 4: The real and imaginary parts of pl 2 (t) versus time 
for the dual avoided crossing model. Same line types as in 
Fig. H (top) P = 50, (bottom) P = 10. 



from the quantum subsystem. The excess coupling term 
in the full mapping QCLE cancels this contribution to 
yield the original QCLE. 

Earlier simulation o 28 ' 35 ' 36 on spin-boson model with bi- 
linear coupling to the bath, along with the present simu- 
lations on two curve crossing models with nonlinear bath 
coupling, indicate that the correction terms are small and 
the approximate PBME yields good agreement with the 
exact quantum results for these models. These results 
suggest that this approximate evolution equation will be 
useful in many applications; however, the validity of the 
PBME must be tested for the specific problem under con- 
sideration and it may not always be a good approxima- 
tion to full QCL dynamics^. The effects of the excess 
coupling term on the time evolution of the observables 
can be estimated by computing the expectation values in 
Eq. ([55)1 and these values provide one way of gauging the 
importance of the deviations from the QCLE. Although 
care must be exercised in the use of the PBME when ap- 
plied to a given system, the ease with which it can be 
simulated and its accuracy for many systems of interest 
make it a powerful simulation scheme that should find 
increased use in future applications. 
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led to an analysis of the PBME, the approximate form of 
the mapping QCLE that retains only the Poisson bracket 
term in the Liouville operator. The PBME is easily simu- 
lated by classical trajectories. However, when compared 
to the full QCLE, it contains an extra term that can 
be interpreted as being proportional to the effect of the 
quantum subsystem on the time evolution of the density 
matrix of the environment along the flow lines generated 
by the classical evolution of the environment in isolation 



Appendix 

In this Appendix we show how exx' (x) in Eq. (|12p and 
<?aa' {%) in Eq. (fTTj) can be computed. 

Calculation of Cx\>(x) 
Starting from its definition we have, 



(2irh) N c X x> (x) = {a\a x ) w 



1 

2h 



dze lp ' z/n (r - -\[q\q\> + pxPx> + i(pxqx' - qxPX')]\r + -) 



= _ / dze lp - z/h 

2h 



l r 
2h 
l r 
2h 
l r 
2h 



dzd(z) 



_d d_ 

dz x dz x > 



• d d . z d . z d 
'd^dz^ + h{r+ 2 )x 'dz- x - hir -2 )x dz^, 
d 



5(z) 



dzy 



*2/* \ f i \ *r Sxx ' l i \ i %( 6xx ' i * V 
n {=P\){tP\') ~ h (— + fr Pxrx '> + h ( Y h Px ' rx > 



r\r\> + i(r\p\> - rypx) +PxPx' ~ fr$\\> 



(47) 
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Similarly for g\y (x) we have 



Calculation of g\y(x) 



UX\>[x)= I dze-^ z l n {m x \r~-){r 



dze- i v*/ h Mri-j)---<Pi(rx 



y) ■■ -M r N ~ -y)<Mn + y)"-0i(rA' + ~^)---M r N + -y) 



1 \ N / 2 2 

h 



(wh) 



dze- lp - z/h \ ( ? AA , e [-3k((-i-^) 2 + (-i + ^) 2 )] . . . ( r 2 _ ^) e [-A((^-^) 2 + (rA + ^) 2 )] 



. . .eh^U^-^^+^+^f) 2 )] + (i _ ^)e [ -*( r?+ ^ )] ■ ■ • (r A - y)e h * ( r ' + ^ )] • • • (r v + y 1 ) 



2 jv+i 

2 N+1 



-x 2 /h 



-x 2 /h 



Sxx'(r 2 x - -(2h - 4pl)) + (1 - 6 X x>)(rx + ipx)(rx> - ipx>) 



rxrx' - i(rxpx> - rx'Px) + Pxpx> - 5xx> 



(48) 



where we have used Eq. together with (f>o(x) — (^-) 1 ^ 4 e 2rX , 4>i(x) = (^) 1//4 y § xe ^ x . In this expression 

recall that z — (z\, Z2, ■ ■ • , zn), P — {pi,P2, ■ ■ ■ ,Pn) and p ■ z — p\Z\ +P2Z2 H +Pn%n, with similar vector notation 

for other unscripted quantities. 
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